if (!requireNamespace("plotly", quietly = TRUE))
    install.packages("plotly")
if (!requireNamespace("kable", quietly = TRUE))
    install.packages("kable")
#> Error in install.packages : error reading from connection
if (!requireNamespace("kableExtra", quietly = TRUE))
    install.packages("kableExtra")

library(comomodels)
library(plotly) # interactive plot
library(knitr) # table layout 
library(kableExtra) # plot and table arrangement

Introduction

This document explains the basics of the SEIaImIsRD model. The model describes how populations of susceptible (\(S\)), exposed (\(E\)), infectious and recovered (\(R\)) individuals evolve over time. The infected population is split into different symptom compartments: asymptomatic (\(I_a\)), mild (\(I_m\)) and severe (\(I_s\)). An ordinary differential equation (ODE) is defined for each of the populations.

The object class SEIaImIsRD lets the user define and simulate such a model, choosing the initial conditions and parameter values.

The SEIaImIsRD Model

Core: the ODE system

The SEIaImIsRD model consists of six ODEs describing the six populations. Susceptible individuals (\(S\)) have not been exposed to the virus. Exposed individuals (\(E\)) have been exposed to the virus and have been infected but are not yet infectious to others; they are in the incubation period. Infectious individuals can spread the virus to others. The total initial population size is normalised to 1. The current model does not include natural death or birth.
The ODEs describing how these six populations evolve over time are \[\frac{\text{d}S}{\text{d}t} = -S (\beta_{I_a}I_a + \beta_{I_m}I_m + \beta_{I_s}I_s) + \omega R,\] where \(\beta\) denotes the rate, at which an infected individual exposes susceptible individuals, and \(\omega\) denotes the rate, at which recovered individuals become susceptible. \[\frac{\text{d}E}{\text{d}t} = S (\beta_{I_a}I_a + \beta_{I_m}I_m + \beta_{I_s}I_s) - \kappa E,\] where \(\kappa\) denotes the rate at which exposed individuals become infected. \[\frac{dI_{a}}{dt} = \eta_{E->I_a} \kappa E - (\gamma_{I_a} + \mu_{I_a})I_{a},\]

\[\frac{dI_{m}}{dt} = \eta_{E->I_m} \kappa E - (\gamma_{I_m} + \mu_{I_a})I_{m},\] \[\frac{dI_{s}}{dt} = \eta_{E->I_s} \kappa E - (\gamma_{I_s} + \mu_{I_a})I_{s},\] where \(\eta_{E->I_a}\), \(\eta_{E->I_m}\), \(\eta_{E->I_s}\) denote the fraction of exposed individuals that develop asymptomatic, mild and symptoms, respectively (in the code the fraction is denoted by ‘p_symptom’). The rate, at which each infected group recovers, is denoted by \(\gamma_I\). The rate of disease-caused mortality of each infected group is denoted by \(\mu_I\).

\[\frac{\text{d}R}{\text{d}t} = - \omega R +\mu_{I_a} * I_{a} + \mu_{I_m} * I_{m} + \mu_{I_s} * I_{s}.\]

The rate at which susceptible individuals are infected depends on the fraction of the population that is susceptible, the fraction of the population that is infectious, and the probability of an infectious individual infecting a susceptible individual, \(\beta\). This is represented in the first terms on the right-hand side of \(\text{d}S/ \text{d}t\) and \(\text{d}I/ \text{d}t\). Exposed individuals become infectious at rate \(\kappa\). The reciprocal value \(1/ \kappa\) is the incubation time of the virus, i.e. how many days after exposure the individual becomes infectious to others. Infectious individuals can either recover or die because of the disease. Infectious individuals in group \(i\in {a,m,s}\) recover in \(1/ \beta_i\) days, thus individuals recover at rate \(\beta_i\). Similarly, infectious individuals die due to the infection after \(1/ \mu_i\) days, making the death rate \(\mu_i\).

In addition to these six main ODEs, the model also keeps track the number of disease-related deaths (\(D\)) over time \[\frac{\text{d}D}{\text{d}t} = \mu_{I_a} * I_{a} + \mu_{I_m} * I_{m} + \mu_{I_s} * I_{s},\]

The system is closed by specifying the initial conditions \[S(0) = S_0,\ E(0) = E_0,\ I(0) = I_0,\ R(0) = R_0,\ C(0) = C_0,\ D(0) = D_0,\] In the implementation in this package, the variables are normalized so \(S(t) + E(t) + I_a(t) + I_m(t) + I_s(t) + R(t) + D(t) \equiv 1\) for any given \(t\).

Steady states of the system described above are computed by setting the right-hand sides of the ODEs equal to zero.

Because we build on this basic model to develop more complex epidemiological models, we solve the ODEs over time instead of just solving for the steady states. The next section shows an example of how to use the SEIaImIsRD model class.

Basic reproduction number

The mathematical details refer to this paper.

We calculate the basic reproduction number (\(\mathscr{R_0}\)) using the next generation matrix approach:

First, define \(x\in (x_1, x_2,...,x_n)^T\) be the number of individuals in each compartment, where the first \(m < n\) compartments contain exposed or infected individuals. In the SEIaImIsRD model, these compartments \(E\), \(I_a\), \(I_m\) and \(I_s\)).

Assume that a disease free equilibrium (DFE) \(x_0\) exists and is stable in the absence of disease, and that the linearized equations for \(x_1, ..., x_m\) at the DFE decouple from the other equations. Consider these equations written in the form \(\frac{\text{d}x_i}{\text{d}t} = \mathscr{F_i}(x) - \mathscr{V_i}(x)\) for \(i \in \{1, 2, ..., m\}\). In this splitting, \(\mathscr{F_i}(x)\) is the rate of appearance of new infections in compartment \(i\), and \(\mathscr{V_i}(x)\) is the rate of other transitions between compartment \(i\) and other infected compartments. It is assumed that \(\mathscr{F_i}=0\) for \(i\in\{m+1, m+2, ..., n\}\).

Now define \(F=\frac{\partial \mathscr{F_i}{(x_0)}}{\partial x_j}\) and \(V=\frac{\partial \mathscr{V_i}{(x_0)}}{\partial x_j}\) for \(i,j \in \{1, 2, ..., m\}\). Matrix \(FV^{-1}\) has \((i,j)\) entry equal to the expected number of secondary infections in compartment \(i\) produced by an infected individual introduced in compartment \(j\). Thus \(FV^{-1}\) is the \(\textit{next generation matrix}\) and

\[\mathscr{R_0} = \rho(FV^{-1}),\] where \(\rho\) denotes the \(\textit{spectral radius}\) of a matrix (i.e. the largest absolute value of its eigenvalues).

Applying the above approach to the SEIaImIsRD model, the DFE matrices \(F\) and \(V\) are: \[ F = \left(\begin{array}{cc} 0 & S_0\beta_{I_a} & S_0\beta_{I_m} & S_0\beta_{I_s}\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 \end{array}\right), V = \left(\begin{array}{cc} \kappa & 0 & 0 & 0\\ -\kappa \eta_{E->I_a} & \gamma_{I_a} + \mu_{I_a} & 0 & 0\\ -\kappa \eta_{E->I_m} & 0 & \gamma_{I_m} + \mu_{I_m} & 0\\ -\kappa \eta_{E->I_s} & 0 & 0 & \gamma_{I_s} + \mu_{I_s} \end{array}\right) \] where \(S_0\) is the initial number of susceptibles, as in assumption the initial conditions \(x_0 = \{S_0, E_0, ...\}\) can be seen as a DFE given that \(S_0, E_0>0\) and \(E_0 \ll S_0\), and \(x_{0i}=0\) for all other compartments except \(S_0\) and \(E_0\) . \(\mathscr{R_0}\) can then be calculated using the equation \(\mathscr{R_0} = \rho(FV^{-1})\).

The SEIaImIsRD class

In this section we will walk through an example of how to use the SEIaImIsRD class.

Initialization of the model

To create a new SEIaImIsRD object, run:

model <- SEIaImIsRD()

Next, we assign values for the fractions of the population. The initial conditions must sum to one.

S <- 0.99
E <- 0.01
I_asymptomatic <- 0
I_mild <- 0
I_severe <- 0
R <- 0
D_cumulative <- 0

Then, we define the transition rates:

beta <- list(i_asymptomatic = 0.80, i_mild = 0.90, i_severe = 1.00)
kappa <- 0.50
omega <- 0.01
p_symptom <- list(i_mild = 0.35, i_severe = 0.05)
gamma <- list(i_asymptomatic = 0.9, i_mild = 0.50, i_severe = 0.05)
mu <- list(i_asymptomatic = 0.005, i_mild = 0.05, i_severe = 0.30)

Define helper function for printing out parameter values

print_table_transmission_parameters <- function(model){
  df.params <- model@transmission_parameters %>%
    data.frame() %>%
    t()
  colnames(df.params) <- "value"
  table.params <- kable(df.params) %>%
    kable_styling(bootstrap_options = "striped", full_width = FALSE, position="left")
  return (table.params)
}

And, we set time for the simulation:

t <- seq(0, 2000, by = 1)

We initialise our model with the parameters defined above:

model <- set_init(object = model, S = S, E = E, I_asymptomatic = I_asymptomatic, I_mild = I_mild,
                    I_severe = I_severe, R = R, D = D_cumulative,
                    beta = beta, kappa = kappa, omega = omega,
                    p_symptom = p_symptom, gamma = gamma, mu = mu)

We use the ode solver to solve the system, plot the relaxation into the steady state and calculate the \(\mathscr{R_0}\) of the system by calling the function R0_SEIaImIsRD. Here we present the simulation result under several different conditions:

1. Including natural loss of immunity

Transmission from \(R\) to \(S\) exists; \(\omega>0\)

par(mfrow=c(1,2))
model <- ode_simulate(model, t)
model <- R0_SEIaImIsRD(model)

# plot
model.plot <- plot_dataframe(model@output, x = "time", y = "fraction", c = "population_group") + 
  ggtitle(paste0("model: R0=", model@R0)) # Show R0 on the title
model.plot <- ggplotly(p=model.plot, dynamicTicks = TRUE, originalData = FALSE)

# get parameters
model.params.table = print_table_transmission_parameters(model)
value
beta.i_asymptomatic 0.800
beta.i_mild 0.900
beta.i_severe 1.000
kappa 0.500
omega 0.010
p_symptom.i_mild 0.350
p_symptom.i_severe 0.050
gamma.i_asymptomatic 0.900
gamma.i_mild 0.500
gamma.i_severe 0.050
mu.i_asymptomatic 0.005
mu.i_mild 0.050
mu.i_severe 0.300

2. No natural loss of immunity

No transmission from \(R\) to \(S\); \(\omega=0\)

model1 <- set_init(object = model, S = S, E = E, I_asymptomatic = I_asymptomatic, I_mild = I_mild,
                    I_severe = I_severe, R = R, D = D_cumulative,
                    beta = beta, kappa = kappa, omega = 0.00,
                    p_symptom = p_symptom, gamma = gamma, mu = mu)

model1 <- ode_simulate(model1, t)
model1 <- R0_SEIaImIsRD(model1)

# plot
model1.plot <- plot_dataframe(model1@output, x = "time", y = "fraction", c = "population_group") +
  ggtitle(paste0("model1: R0=", model1@R0))
model1.plot <- ggplotly(p=model1.plot, dynamicTicks = TRUE, originalData = FALSE)

# get parameters
model1.params.table = print_table_transmission_parameters(model1)
value
beta.i_asymptomatic 0.800
beta.i_mild 0.900
beta.i_severe 1.000
kappa 0.500
omega 0.000
p_symptom.i_mild 0.350
p_symptom.i_severe 0.050
gamma.i_asymptomatic 0.900
gamma.i_mild 0.500
gamma.i_severe 0.050
mu.i_asymptomatic 0.005
mu.i_mild 0.050
mu.i_severe 0.300

We can see whether including waning immunity may produce different epidemic behaviors under the same \(\mathscr{R_0}\).

3. Altering the value(s) of \(\beta\)

\(\beta\): the rate at which an individual from any infected group \(I_i\) infects \(S\)

beta <- list(i_asymptomatic = 0.50, i_mild = 0.90, i_severe = 1.00)
model2 <- set_init(object = model, S = S, E = E, I_asymptomatic = I_asymptomatic, I_mild = I_mild,
                    I_severe = I_severe, R = R, D = D_cumulative,
                    beta = beta, kappa = kappa, omega = 0.00,
                    p_symptom = p_symptom, gamma = gamma, mu = mu)

model2 <- ode_simulate(model2, t)
model2 <- R0_SEIaImIsRD(model2)

# plot
model2.plot <- plot_dataframe(model2@output, x = "time", y = "fraction", c = "population_group") +
  ggtitle(paste0("model2: R0=", model1@R0))
model2.plot <- ggplotly(p=model2.plot, dynamicTicks = TRUE, originalData = FALSE)

# get parameters
model2.params.table = print_table_transmission_parameters(model2)
value
beta.i_asymptomatic 0.500
beta.i_mild 0.900
beta.i_severe 1.000
kappa 0.500
omega 0.000
p_symptom.i_mild 0.350
p_symptom.i_severe 0.050
gamma.i_asymptomatic 0.900
gamma.i_mild 0.500
gamma.i_severe 0.050
mu.i_asymptomatic 0.005
mu.i_mild 0.050
mu.i_severe 0.300